# Model one 
model{	
	for(i in 1:N){	# loop over observations 

	   # Utility Decision 1
	   U11[i] <- beta[1,1] + delta[1] * (gamma[1] * Up1[i] + (1-gamma[1]) * Ucp[i])  + delta[6] *lr1[i]  + delta[3] * PID1[i]  + X[i,] %*% psi[1,1,1:6] 
	   U21[i] <- beta[1,2] + delta[2] * (gamma[2] * Up1[i] + (1-gamma[2]) * Ucp[i])  + delta[7] *lr1[i]  + delta[4] * PID1[i]  + delta[5] * a11[i] + X[i,] %*% psi[1,2,1:6] 

	   # Utility Decision 2
	   U12[i] <- beta[2,1] + delta[1] * (gamma[1] * Up2[i] + (1-gamma[1]) * Ucp[i])  + delta[6] *lr2[i] + delta[3] * PID2[i] + X[i,] %*% psi[2,1,1:6] 
	   U22[i] <- beta[2,2] + delta[2] * (gamma[2] * Up2[i] + (1-gamma[2]) * Ucp[i])  + delta[7] *lr2[i]  + delta[4] * PID2[i] + delta[5] * a12[i] + X[i,] %*% psi[2,2,1:6] 


	   # Utility 3
	   U13[i] <- 0
	   U23[i] <- 0

	    # Normalize

	    P1[i,1] <- exp(U11[i]) / (exp(U11[i]) + exp(U12[i]) + exp(U13[i]))  
	    P1[i,2] <- exp(U12[i]) / (exp(U11[i]) + exp(U12[i]) + exp(U13[i]))  
	    P1[i,3] <- exp(U13[i]) / (exp(U11[i]) + exp(U12[i]) + exp(U13[i]))  

	    P2[i,1] <- exp(U21[i]) / (exp(U21[i]) + exp(U22[i]) + exp(U23[i]))  
	    P2[i,2] <- exp(U22[i]) / (exp(U21[i]) + exp(U22[i]) + exp(U23[i]))  
	    P2[i,3] <- exp(U23[i]) / (exp(U21[i]) + exp(U22[i]) + exp(U23[i]))  


	   # Probailities
	   P[i,1] <-  P1[i,1] * P2[i,1]
	   P[i,2] <-  P1[i,1] * P2[i,2]
	   P[i,3] <-  P1[i,1] * P2[i,3]

	   P[i,4] <-  P1[i,2] * P2[i,1]
	   P[i,5] <-  P1[i,2] * P2[i,2]
	   P[i,6] <-  P1[i,2] * P2[i,3]

	   P[i,7] <-  P1[i,3] * P2[i,1]
	   P[i,8] <-  P1[i,3] * P2[i,2]
	   P[i,9] <-  P1[i,3] * P2[i,3]


	   # Distribution
	   Y[i] ~ dcat(P[i,1:9])

       }

       	 	  
       gamma[1] ~ dunif(0,1)
       gamma[2] ~ dunif(0,1)
	


	# Other Prior
	delta[1] ~ dnorm(0,.001)
	delta[2] ~ dnorm(0,.001)
	delta[3] ~ dnorm(0,.001)
	delta[4] ~ dnorm(0,.001)
	delta[5] ~ dnorm(0,.001)
	delta[6] ~ dnorm(0,.001)
	delta[7] ~ dnorm(0,.001)

	beta[1,1] ~ dnorm(0,0.001)
	beta[1,2] ~ dnorm(0,0.001)
	beta[2,1] ~ dnorm(0,0.001)
	beta[2,2] ~ dnorm(0,0.001)
	
	# psi
	for(k in 1:6){
	psi[1,1,k]  ~ dnorm(0,.001)
	psi[1,2,k]  ~ dnorm(0,.001)
	psi[2,1,k]  ~ dnorm(0,.001)
	psi[2,2,k]  ~ dnorm(0,.001)		
	}

	# Sample Difference
	diff <- gamma[1] - gamma[2]

}
